home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 22 / Cream of the Crop 22.iso / math / ast53src.zip / CHARTS3.C < prev    next >
C/C++ Source or Header  |  1996-09-29  |  26KB  |  715 lines
  1. /*
  2. ** Astrolog (Version 5.30) File: charts3.c
  3. **
  4. ** IMPORTANT NOTICE: The graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1996 by Walter D. Pullen
  6. ** (Astara@msn.com, http://www.magitech.com/~cruiser1/astrolog.htm).
  7. ** Permission is granted to freely use and distribute these routines
  8. ** provided one doesn't sell, restrict, or profit from them in any way.
  9. ** Modification is allowed provided these notices remain with any
  10. ** altered or edited versions of the program.
  11. **
  12. ** The main planetary calculation routines used in this program have
  13. ** been Copyrighted and the core of this program is basically a
  14. ** conversion to C of the routines created by James Neely as listed in
  15. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  16. ** available from Matrix Software. The copyright gives us permission to
  17. ** use the routines for personal use but not to sell them or profit from
  18. ** them in any way.
  19. **
  20. ** The PostScript code within the core graphics routines are programmed
  21. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  22. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  23. **
  24. ** The extended accurate ephemeris databases and formulas are from the
  25. ** calculation routines in the program "Placalc" and are programmed and
  26. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  27. ** (alois@azur.ch). The use of that source code is subject to
  28. ** regulations made by Astrodienst Zurich, and the code is not in the
  29. ** public domain. This copyright notice must not be changed or removed
  30. ** by any user of this program.
  31. **
  32. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  33. ** X Window graphics initially programmed 10/23-29/1991.
  34. ** PostScript graphics initially programmed 11/29-30/1992.
  35. ** Last code change made 9/22/1996.
  36. */
  37.  
  38. #include "astrolog.h"
  39.  
  40.  
  41. /*
  42. ******************************************************************************
  43. ** Multiple Chart Scanning Routines.
  44. ******************************************************************************
  45. */
  46.  
  47. /* Search through a day, and print out the times of exact aspects among the  */
  48. /* planets during that day, as specified with the -d switch, as well as the  */
  49. /* times when a planet changes sign or direction. To do this, we cast charts */
  50. /* for the beginning and end of the day, or a part of a day, and do a linear */
  51. /* equation check to see if anything exciting happens during the interval.   */
  52. /* (This is probably the single most complicated procedure in the program.)  */
  53.  
  54. void ChartInDaySearch(fProg)
  55. bool fProg;
  56. {
  57.   char sz[cchSzDef];
  58.   int source[MAXINDAY], aspect[MAXINDAY], dest[MAXINDAY],
  59.     sign1[MAXINDAY], sign2[MAXINDAY], D1, D2, counttotal = 0, occurcount,
  60.     division, div, fYear, yea0, yea1, yea2, i, j, k, l, s1, s2;
  61.   real time[MAXINDAY], divsiz, d1, d2, e1, e2, f1, f2, g;
  62.   CI ciT;
  63.  
  64.   /* If parameter 'fProg' is set, look for changes in a progressed chart. */
  65.  
  66.   ciT = ciTran;
  67.   fYear = us.fInDayMonth && (MonT == 0);
  68.   division = (fYear || fProg) ? 1 : us.nDivision;
  69.   divsiz = 24.0 / (real)division*60.0;
  70.  
  71.   /* If -dY in effect, then search through a range of years. */
  72.  
  73.   yea1 = fProg ? YeaT : Yea;
  74.   yea2 = fYear ? (yea1 + us.nEphemYears - 1) : yea1;
  75.   for (yea0 = yea1; yea0 <= yea2; yea0++) {
  76.  
  77.   /* If -dm in effect, then search through the whole month, day by day. */
  78.  
  79.   if (us.fInDayMonth) {
  80.     D1 = 1;
  81.     if (fYear) {
  82.       MonT = 1; D2 = DayInYearHi(yea0);
  83.     } else
  84.       D2 = DayInMonth(fProg ? MonT : Mon, yea0);
  85.   } else
  86.     D1 = D2 = Day;
  87.  
  88.   /* Start searching the day or days in question for exciting stuff. */
  89.  
  90.   for (DayT = D1; DayT <= D2; DayT = AddDay(Mon, DayT, yea0, 1)) {
  91.     occurcount = 0;
  92.  
  93.     /* Cast chart for beginning of day and store it for future use. */
  94.  
  95.     SetCI(ciCore, fYear ? MonT : Mon, DayT, yea0, 0.0, Dst, Zon, Lon, Lat);
  96.     if (us.fProgress = fProg) {
  97.       is.JDp = MdytszToJulian(MonT, DD, yea0, 0.0, Dst, Zon);
  98.       ciCore = ciMain;
  99.     }
  100.     CastChart(fTrue);
  101.     cp2 = cp0;
  102.  
  103.     /* Now divide the day into segments and search each segment in turn. */
  104.     /* More segments is slower, but has slightly better time accuracy.   */
  105.  
  106.     for (div = 1; div <= division; div++) {
  107.  
  108.       /* Cast the chart for the ending time of the present segment. The   */
  109.       /* beginning time chart is copied from the previous end time chart. */
  110.  
  111.       SetCI(ciCore, fYear ? MonT : Mon, DayT, yea0,
  112.         DegToDec(24.0*(real)div/(real)division), Dst, Zon, Lon, Lat);
  113.       if (fProg) {
  114.         is.JDp = MdytszToJulian(MonT, DD+1, yea0, 0.0, Dst, Zon);
  115.         ciCore = ciMain;
  116.       }
  117.       CastChart(fTrue);
  118.       cp1 = cp2; cp2 = cp0;
  119.  
  120.       /* Now search through the present segment for anything exciting. */
  121.  
  122.       for (i = 0; i <= cObj; i++) if (!FIgnore(i) && (fProg || FThing(i))) {
  123.         s1 = SFromZ(cp1.obj[i])-1;
  124.         s2 = SFromZ(cp2.obj[i])-1;
  125.  
  126.         /* Does the current planet change into the next or previous sign? */
  127.  
  128.         if (s1 != s2 && !us.fIgnoreSign) {
  129.           source[occurcount] = i;
  130.           aspect[occurcount] = aSig;
  131.           dest[occurcount] = s2+1;
  132.           time[occurcount] = MinDistance(cp1.obj[i],
  133.             (real)(cp1.dir[i] >= 0.0 ? s2 : s1) * 30.0) /
  134.             MinDistance(cp1.obj[i], cp2.obj[i])*divsiz + (real)(div-1)*divsiz;
  135.           sign1[occurcount] = sign2[occurcount] = s1+1;
  136.           occurcount++;
  137.         }
  138.  
  139.         /* Does the current planet go retrograde or direct? */
  140.  
  141.         if ((cp1.dir[i] < 0.0) != (cp2.dir[i] < 0.0) && !us.fIgnoreDir) {
  142.           source[occurcount] = i;
  143.           aspect[occurcount] = aDir;
  144.           dest[occurcount] = cp2.dir[i] < 0.0;
  145.           time[occurcount] = RAbs(cp1.dir[i])/(RAbs(cp1.dir[i])+
  146.             RAbs(cp2.dir[i]))*divsiz + (real)(div-1)*divsiz;
  147.           sign1[occurcount] = sign2[occurcount] = s1+1;
  148.           occurcount++;
  149.         }
  150.  
  151.         /* Now search for anything making an aspect to the current planet. */
  152.  
  153.         for (j = i+1; j <= cObj; j++) if (!FIgnore(j) && (fProg || FThing(j)))
  154.           for (k = 1; k <= us.nAsp; k++) if (FAcceptAspect(i, k, j)) {
  155.             d1 = cp1.obj[i]; d2 = cp2.obj[i];
  156.             e1 = cp1.obj[j]; e2 = cp2.obj[j];
  157.             if (MinDistance(d1, d2) < MinDistance(e1, e2)) {
  158.               SwapR(&d1, &e1);
  159.               SwapR(&d2, &e2);
  160.             }
  161.  
  162.             /* We are searching each aspect in turn. Let's subtract the  */
  163.             /* size of the aspect from the angular difference, so we can */
  164.             /* then treat it like a conjunction.                         */
  165.  
  166.             if (MinDistance(e1, Mod(d1-rAspAngle[k])) <
  167.                 MinDistance(e2, Mod(d2+rAspAngle[k]))) {
  168.               e1 = Mod(e1+rAspAngle[k]);
  169.               e2 = Mod(e2+rAspAngle[k]);
  170.             } else {
  171.               e1 = Mod(e1-rAspAngle[k]);
  172.               e2 = Mod(e2-rAspAngle[k]);
  173.             }
  174.  
  175.             /* Check to see if the aspect actually occurs during our    */
  176.             /* segment, making sure we take into account if one or both */
  177.             /* planets are retrograde or if they cross the Aries point. */
  178.  
  179.             f1 = e1-d1;
  180.             if (RAbs(f1) > rDegHalf)
  181.               f1 -= RSgn(f1)*rDegMax;
  182.             f2 = e2-d2;
  183.             if (RAbs(f2) > rDegHalf)
  184.               f2 -= RSgn(f2)*rDegMax;
  185.             if (MinDistance(Midpoint(d1, d2), Midpoint(e1, e2)) < rDegQuad &&
  186.               RSgn(f1) != RSgn(f2)) {
  187.               source[occurcount] = i;
  188.               aspect[occurcount] = k;
  189.               dest[occurcount] = j;
  190.  
  191.               /* Horray! The aspect occurs sometime during the interval.   */
  192.               /* Now we just have to solve an equation in two variables to */
  193.               /* find out where the "lines" cross, i.e. the aspect's time. */
  194.  
  195.               f1 = d2-d1;
  196.               if (RAbs(f1) > rDegHalf)
  197.                 f1 -= RSgn(f1)*rDegMax;
  198.               f2 = e2-e1;
  199.               if (RAbs(f2) > rDegHalf)
  200.                 f2 -= RSgn(f2)*rDegMax;
  201.               g = (RAbs(d1-e1) > rDegHalf ?
  202.                 (d1-e1)-RSgn(d1-e1)*rDegMax : d1-e1)/(f2-f1);
  203.               time[occurcount] = g*divsiz + (real)(div-1)*divsiz;
  204.               sign1[occurcount] = (int)(Mod(cp1.obj[i]+
  205.                 RSgn(cp2.obj[i]-cp1.obj[i])*
  206.                 (RAbs(cp2.obj[i]-cp1.obj[i]) > rDegHalf ? -1 : 1)*
  207.                 RAbs(g)*MinDistance(cp1.obj[i], cp2.obj[i]))/30.0)+1;
  208.               sign2[occurcount] = (int)(Mod(cp1.obj[j]+
  209.                 RSgn(cp2.obj[j]-cp1.obj[j])*
  210.                 (RAbs(cp2.obj[j]-cp1.obj[j]) > rDegHalf ? -1 : 1)*
  211.                 RAbs(g)*MinDistance(cp1.obj[j], cp2.obj[j]))/30.0)+1;
  212.               occurcount++;
  213.             }
  214.           }
  215.       }
  216.     }
  217.  
  218.     /* After all the aspects, etc, in the day have been located, sort   */
  219.     /* them by time at which they occur, so we can print them in order. */
  220.  
  221.     for (i = 1; i < occurcount; i++) {
  222.       j = i-1;
  223.       while (j >= 0 && time[j] > time[j+1]) {
  224.         SwapN(source[j], source[j+1]);
  225.         SwapN(aspect[j], aspect[j+1]);
  226.         SwapN(dest[j], dest[j+1]);
  227.         SwapR(&time[j], &time[j+1]);
  228.         SwapN(sign1[j], sign1[j+1]); SwapN(sign2[j], sign2[j+1]);
  229.         j--;
  230.       }
  231.     }
  232.  
  233.     /* Finally, loop through and display each aspect and when it occurs. */
  234.  
  235.     for (i = 0; i < occurcount; i++) {
  236.       s1 = (int)time[i]/60;
  237.       s2 = (int)time[i]-s1*60;
  238.       j = DayT;
  239.       if (fYear || fProg) {
  240.         l = MonT;
  241.         while (j > (k = DayInMonth(l, yea0))) {
  242.           j -= k;
  243.           l++;
  244.         }
  245.       }
  246.       SetCI(ciSave, fYear || fProg ? l : Mon, j, yea0,
  247.         DegToDec(time[i] / 60.0), Dst, Zon, Lon, Lat);
  248.       k = DayOfWeek(fYear || fProg ? l : Mon, j, yea0);
  249.       AnsiColor(kRainbowA[k + 1]);
  250.       sprintf(sz, "(%c%c%c) ", chDay3(k)); PrintSz(sz);
  251.       AnsiColor(kDefault);
  252.       sprintf(sz, "%s %s ",
  253.         SzDate(fYear || fProg ? l : Mon, j, yea0, fFalse),
  254.         SzTime(s1, s2)); PrintSz(sz);
  255.       PrintAspect(source[i], sign1[i],
  256.         (int)RSgn(cp1.dir[source[i]])+(int)RSgn(cp2.dir[source[i]]),
  257.         aspect[i], dest[i], sign2[i],
  258.         (int)RSgn(cp1.dir[dest[i]])+(int)RSgn(cp2.dir[dest[i]]),
  259.         (char)(fProg ? 'e' : 'd'));
  260.       PrintInDay(source[i], aspect[i], dest[i]);
  261.     }
  262.     counttotal += occurcount;
  263.   }
  264.   }
  265.   if (counttotal == 0)
  266.     PrintSz("No transit events found.\n");
  267.  
  268.   /* Recompute original chart placements as we've overwritten them. */
  269.  
  270.   ciCore = ciMain; ciTran = ciT;
  271.   CastChart(fTrue);
  272. }
  273.  
  274.  
  275. /* Search through a month, year, or years, and print out the times of exact */
  276. /* transits where planets in the time frame make aspect to the planets in   */
  277. /* some other chart, as specified with the -t switch. To do this, we cast   */
  278. /* charts for the start and end of each month, or within a month, and do an */
  279. /* equation check for aspects to the other base chart during the interval.  */
  280.  
  281. void ChartTransitSearch(fProg)
  282. bool fProg;
  283. {
  284.   real planet3[objMax], house3[cSign+1], ret3[objMax], time[MAXINDAY];
  285.   char sz[cchSzDef];
  286.   int source[MAXINDAY], aspect[MAXINDAY], dest[MAXINDAY], sign[MAXINDAY],
  287.     isret[MAXINDAY], M1, M2, Y1, Y2, counttotal = 0, occurcount, division,
  288.     div, nAsp, fCusp, i, j, k, s1, s2, s3;
  289.   real divsiz, daysiz, d, e1, e2, f1, f2;
  290.   CI ciT;
  291.  
  292.   /* Save away natal chart and initialize things. */
  293.  
  294.   ciT = ciTran;
  295.   for (i = 1; i <= cSign; i++)
  296.     house3[i] = chouse[i];
  297.   for (i = 0; i <= cObj; i++) {
  298.     planet3[i] = planet[i];
  299.     ret3[i] = ret[i];
  300.   }
  301.   if (fProg)
  302.     fCusp = fFalse;
  303.   else {
  304.     fCusp = fTrue;
  305.     for (i = cuspLo; i <= cuspHi; i++)
  306.       fCusp &= ignore2[i];
  307.   }
  308.   division = us.nDivision;
  309.   if (!fProg && !fCusp)
  310.     division = Max(division, 96);
  311.   nAsp = is.fReturn ? aCon : us.nAsp;
  312.  
  313.   /* Hacks: Searching month number zero means to search the whole year    */
  314.   /* instead, month by month. Searching a negative month means to search  */
  315.   /* multiple years, with the span of the year stored in the "day" field. */
  316.  
  317.   Y1 = Y2 = YeaT;
  318.   M1 = M2 = MonT;
  319.   if (MonT < 1) {
  320.     M1 = 1; M2 = 12;
  321.     if (MonT < 0) {
  322.       if (DayT < 1) {
  323.         Y1 = YeaT + DayT + 1; Y2 = YeaT;
  324.       } else {
  325.         Y1 = YeaT; Y2 = YeaT + DayT - 1;
  326.       }
  327.     }
  328.   }
  329.  
  330.   /* Start searching the year or years in question for any transits. */
  331.  
  332.   for (YeaT = Y1; YeaT <= Y2; YeaT++)
  333.  
  334.   /* Start searching the month or months in question for any transits. */
  335.  
  336.   for (MonT = M1; MonT <= M2; MonT++) {
  337.     daysiz = (real)DayInMonth(MonT, YeaT)*24.0*60.0;
  338.     divsiz = daysiz / (real)division;
  339.  
  340.     /* Cast chart for beginning of month and store it for future use. */
  341.  
  342.     SetCI(ciCore, MonT, 1, YeaT, 0.0, DstT, ZonT, LonT, LatT);
  343.     if (us.fProgress = fProg) {
  344.       is.JDp = MdytszToJulian(MM, DD, YY, 0.0, DstT, ZonT);
  345.       ciCore = ciMain;
  346.     }
  347.     for (i = 0; i <= oNorm; i++)
  348.       SwapN(ignore[i], ignore2[i]);
  349.     CastChart(fTrue);
  350.     for (i = 0; i <= oNorm; i++)
  351.       SwapN(ignore[i], ignore2[i]);
  352.     cp2 = cp0;
  353.  
  354.     /* Divide our month into segments and then search each segment in turn. */
  355.  
  356.     for (div = 1; div <= division; div++) {
  357.       occurcount = 0;
  358.  
  359.       /* Cast the chart for the ending time of the present segment, and */
  360.       /* copy the start time chart from the previous end time chart.    */
  361.  
  362.       d = 1.0 + (daysiz/24.0/60.0)*(real)div/(real)division;
  363.       SetCI(ciCore, MonT, (int)d, YeaT,
  364.         DegToDec(RFract(d)*24.0), DstT, ZonT, LonT, LatT);
  365.       if (fProg) {
  366.         is.JDp = MdytszToJulian(MM, DD, YY, 0.0, DstT, ZonT);
  367.         ciCore = ciMain;
  368.       }
  369.       for (i = 0; i <= oNorm; i++)
  370.         SwapN(ignore[i], ignore2[i]);
  371.       CastChart(fTrue);
  372.       for (i = 0; i <= oNorm; i++)
  373.         SwapN(ignore[i], ignore2[i]);
  374.       cp1 = cp2; cp2 = cp0;
  375.  
  376.       /* Now search through the present segment for any transits. Note that */
  377.       /* stars can be transited, but they can't make transits themselves.   */
  378.  
  379.       for (i = 0; i <= cObj; i++) if (!FIgnore(i)) {
  380.         for (j = 0; j <= oNorm; j++) {
  381.           if ((is.fReturn ? i != j : FIgnore2(j)) || (fCusp && !FThing(j)))
  382.             continue;
  383.  
  384.           /* Between each pair of planets, check if they make any aspects. */
  385.  
  386.           for (k = 1; k <= nAsp; k++) if (FAcceptAspect(i, k, j)) {
  387.             d = planet3[i]; e1 = cp1.obj[j]; e2 = cp2.obj[j];
  388.             if (MinDistance(e1, Mod(d-rAspAngle[k])) <
  389.                 MinDistance(e2, Mod(d+rAspAngle[k]))) {
  390.               e1 = Mod(e1+rAspAngle[k]);
  391.               e2 = Mod(e2+rAspAngle[k]);
  392.             } else {
  393.               e1 = Mod(e1-rAspAngle[k]);
  394.               e2 = Mod(e2-rAspAngle[k]);
  395.             }
  396.  
  397.             /* Check to see if the present aspect actually occurs during the */
  398.             /* segment, making sure we check any Aries point crossings.      */
  399.  
  400.             f1 = e1-d;
  401.             if (RAbs(f1) > rDegHalf)
  402.               f1 -= RSgn(f1)*rDegMax;
  403.             f2 = e2-d;
  404.             if (RAbs(f2) > rDegHalf)
  405.               f2 -= RSgn(f2)*rDegMax;
  406.             if (MinDistance(d, Midpoint(e1, e2)) < rDegQuad &&
  407.               RSgn(f1) != RSgn(f2) && occurcount < MAXINDAY) {
  408.  
  409.               /* Ok, we have found a transit. Now determine the time */
  410.               /* and save this transit in our list to be printed.    */
  411.  
  412.               source[occurcount] = j;
  413.               aspect[occurcount] = k;
  414.               dest[occurcount] = i;
  415.               time[occurcount] = RAbs(f1)/(RAbs(f1)+RAbs(f2))*divsiz +
  416.                 (real)(div-1)*divsiz;
  417.               sign[occurcount] = (int)(Mod(
  418.                 MinDistance(cp1.obj[j], Mod(d-rAspAngle[k])) <
  419.                 MinDistance(cp2.obj[j], Mod(d+rAspAngle[k])) ?
  420.                 d-rAspAngle[k] : d+rAspAngle[k])/30.0)+1;
  421.               isret[occurcount] = (int)RSgn(cp1.dir[j]) +
  422.                 (int)RSgn(cp2.dir[j]);
  423.               occurcount++;
  424.             }
  425.           }
  426.         }
  427.       }
  428.  
  429.       /* After all transits located, sort them by time at which they occur. */
  430.  
  431.       for (i = 1; i < occurcount; i++) {
  432.         j = i-1;
  433.         while (j >= 0 && time[j] > time[j+1]) {
  434.           SwapN(source[j], source[j+1]);
  435.           SwapN(aspect[j], aspect[j+1]);
  436.           SwapN(dest[j], dest[j+1]);
  437.           SwapR(&time[j], &time[j+1]);
  438.           SwapN(sign[j], sign[j+1]);
  439.           SwapN(isret[j], isret[j+1]);
  440.           j--;
  441.         }
  442.       }
  443.  
  444.       /* Now loop through list and display all the transits. */
  445.  
  446.       for (i = 0; i < occurcount; i++) {
  447.         s1 = (_int)time[i]/24/60;
  448.         s3 = (_int)time[i]-s1*24*60;
  449.         s2 = s3/60;
  450.         s3 = s3-s2*60;
  451.         SetCI(ciSave, MonT, s1+1, YeaT, DegToDec((real)
  452.           ((_int)time[i]-s1*24*60) / 60.0), DstT, ZonT, LonT, LatT);
  453.         sprintf(sz, "%s %s ",
  454.           SzDate(MonT, s1+1, YeaT, fFalse), SzTime(s2, s3)); PrintSz(sz);
  455.         PrintAspect(source[i], sign[i], isret[i], aspect[i],
  456.           dest[i], SFromZ(planet3[dest[i]]), (int)RSgn(ret3[dest[i]]),
  457.           (char)(fProg ? 'u' : 't'));
  458.  
  459.         /* Check for a Solar, Lunar, or any other return. */
  460.  
  461.         if (aspect[i] == aCon && source[i] == dest[i]) {
  462.           AnsiColor(kWhite);
  463.           sprintf(sz, " (%s Return)", source[i] == oSun ? "Solar" :
  464.             (source[i] == oMoo ? "Lunar" : szObjName[source[i]]));
  465.           PrintSz(sz);
  466.         }
  467.         PrintL();
  468. #ifdef INTERPRET
  469.         if (us.fInterpret)
  470.           InterpretTransit(source[i], aspect[i], dest[i]);
  471. #endif
  472.         AnsiColor(kDefault);
  473.       }
  474.       counttotal += occurcount;
  475.     }
  476.   }
  477.   if (counttotal == 0)
  478.     PrintSz("No transits found.\n");
  479.  
  480.   /* Recompute original chart placements as we've overwritten them. */
  481.  
  482.   ciCore = ciMain; ciTran = ciT;
  483.   us.fProgress = fFalse;
  484.   CastChart(fTrue);
  485. }
  486.  
  487.  
  488. /* Display a list of planetary rising times relative to the local horizon */
  489. /* for the day indicated in the chart information, as specified with the  */
  490. /* -Zd switch. For the day, the time each planet rises (transits horizon  */
  491. /* in East half of sky), sets (transits horizon in West), reaches its     */
  492. /* zenith point (transits meridian in South half of sky), and nadirs      */
  493. /* transits meridian in North), is displayed.                             */
  494.  
  495. void ChartInDayHorizon()
  496. {
  497.   char sz[cchSzDef];
  498.   int source[MAXINDAY], type[MAXINDAY], sign[MAXINDAY],
  499.     fRet[MAXINDAY], occurcount, division, div, s1, s2, i, j, fT;
  500.   real time[MAXINDAY], rgalt1[objMax], rgalt2[objMax], azialt[MAXINDAY],
  501.     azi1, azi2, alt1, alt2, lon, lat, mc1, mc2, xA, yA, xV, yV, d, k;
  502.   CI ciT;
  503.  
  504.   fT = us.fSidereal; us.fSidereal = fFalse;
  505.   lon = RFromD(Mod(Lon)); lat = RFromD(Lat);
  506.   division = us.nDivision * 4;
  507.   occurcount = 0;
  508.  
  509.   ciT = ciTwin; ciCore = ciMain; ciCore.tim = 0.0;
  510.   CastChart(fTrue);
  511.   mc2 = RFromD(planet[oMC]); k = RFromD(planetalt[oMC]);
  512.   EclToEqu(&mc2, &k);
  513.   cp2 = cp0;
  514.   for (i = 1; i <= cObj; i++) {
  515.     rgalt2[i] = planetalt[i];
  516.   }
  517.  
  518.   /* Loop through the day, dividing it into a certain number of segments. */
  519.   /* For each segment we get the planet positions at its endpoints.       */
  520.  
  521.   for (div = 1; div <= division; div++) {
  522.     ciCore = ciMain; ciCore.tim = DegToDec(24.0*(real)div/(real)division);
  523.     CastChart(fTrue);
  524.     mc1 = mc2;
  525.     mc2 = RFromD(planet[oMC]); k = RFromD(planetalt[oMC]);
  526.     EclToEqu(&mc2, &k);
  527.     cp1 = cp2; cp2 = cp0;
  528.     for (i = 1; i <= cObj; i++) {
  529.       rgalt1[i] = rgalt2[i]; rgalt2[i] = planetalt[i];
  530.     }
  531.  
  532.     /* For our segment, check to see if each planet during it rises, sets, */
  533.     /* reaches its zenith, or reaches its nadir.                           */
  534.  
  535.     for (i = 1; i <= cObj; i++) if (!ignore[i] && FThing(i)) {
  536.       EclToHorizon(&azi1, &alt1, cp1.obj[i], rgalt1[i], lon, lat, mc1);
  537.       EclToHorizon(&azi2, &alt2, cp2.obj[i], rgalt2[i], lon, lat, mc2);
  538.       j = 0;
  539.  
  540.       /* Check for transits to the horizon. */
  541.       if ((alt1 > 0.0) != (alt2 > 0.0)) {
  542.         d = RAbs(alt1)/(RAbs(alt1)+RAbs(alt2));
  543.         k = Mod(azi1 + d*MinDifference(azi1, azi2));
  544.         j = 1 + 2*(MinDistance(k, rDegHalf) < rDegQuad);
  545.  
  546.       /* Check for transits to the meridian. */
  547.       } else if (RSgn(MinDifference(azi1, rDegQuad)) !=
  548.         RSgn(MinDifference(azi2, rDegQuad))) {
  549.         j = 2 + 2*(MinDistance(azi1, rDegQuad) < rDegQuad);
  550.         d = RAbs(azi1 - (j > 2 ? rDegQuad : 270.0))/MinDistance(azi1, azi2);
  551.         k = alt1 + d*(alt2-alt1);
  552.       }
  553.       if (j && !ignorez[j-1] && occurcount < MAXINDAY) {
  554.         source[occurcount] = i;
  555.         type[occurcount] = j;
  556.         time[occurcount] = 24.0*((real)(div-1)+d)/(real)division*60.0;
  557.         sign[occurcount] = (int)Mod(cp1.obj[i] +
  558.           d*MinDifference(cp1.obj[i], cp2.obj[i]))/30 + 1;
  559.         fRet[occurcount] = (int)RSgn(cp1.dir[i]) + (int)RSgn(cp2.dir[i]);
  560.         azialt[occurcount] = k;
  561.         ciSave = ciMain;
  562.         ciSave.tim = DegToDec(time[occurcount] / 60.0);
  563.         occurcount++;
  564.       }
  565.     }
  566.   }
  567.  
  568.   /* Sort each event in order of time when it happens during the day. */
  569.  
  570.   for (i = 1; i < occurcount; i++) {
  571.     j = i-1;
  572.     while (j >= 0 && time[j] > time[j+1]) {
  573.       SwapN(source[j], source[j+1]);
  574.       SwapN(type[j], type[j+1]);
  575.       SwapR(&time[j], &time[j+1]);
  576.       SwapN(sign[j], sign[j+1]);
  577.       SwapN(fRet[j], fRet[j+1]);
  578.       SwapR(&azialt[j], &azialt[j+1]);
  579.       j--;
  580.     }
  581.   }
  582.  
  583.   /* Finally display the list showing each event and when it occurs. */
  584.  
  585.   for (i = 0; i < occurcount; i++) {
  586.     ciSave = ciMain;
  587.     ciSave.tim = DegToDec(time[i] / 60.0);
  588.     j = DayOfWeek(Mon, Day, Yea);
  589.     AnsiColor(kRainbowA[j + 1]);
  590.     sprintf(sz, "(%c%c%c) ", chDay3(j)); PrintSz(sz);
  591.     AnsiColor(kDefault);
  592.     s1 = (int)time[i]/60;
  593.     s2 = (int)time[i]-s1*60;
  594.     sprintf(sz, "%s %s ", SzDate(Mon, Day, Yea, fFalse), SzTime(s1, s2));
  595.     PrintSz(sz);
  596.     AnsiColor(kObjA[source[i]]);
  597.     sprintf(sz, "%7.7s ", szObjName[source[i]]); PrintSz(sz);
  598.     AnsiColor(kSignA(sign[i]));
  599.     sprintf(sz, "%c%c%c%c%c ",
  600.       fRet[i] > 0 ? '(' : (fRet[i] < 0 ? '[' : '<'), chSig3(sign[i]),
  601.       fRet[i] > 0 ? ')' : (fRet[i] < 0 ? ']' : '>')); PrintSz(sz);
  602.     AnsiColor(kElemA[type[i]-1]);
  603.     if (type[i] == 1)
  604.       PrintSz("rises  ");
  605.     else if (type[i] == 2)
  606.       PrintSz("zeniths");
  607.     else if (type[i] == 3)
  608.       PrintSz("sets   ");
  609.     else
  610.       PrintSz("nadirs ");
  611.     AnsiColor(kDefault);
  612.     PrintSz(" at ");
  613.     if (type[i] & 1) {
  614.       j = (int)(RFract(azialt[i])*60.0);
  615.       sprintf(sz, "%3d%c%02d'", (int)azialt[i], chDeg1, j); PrintSz(sz);
  616.  
  617.       /* For rising and setting events, we'll also display a direction  */
  618.       /* vector to make the 360 degree azimuth value thought of easier. */
  619.  
  620.       xA = RCosD(azialt[i]); yA = RSinD(azialt[i]);
  621.       if (RAbs(xA) < RAbs(yA)) {
  622.         xV = RAbs(xA / yA); yV = 1.0;
  623.       } else {
  624.         yV = RAbs(yA / xA); xV = 1.0;
  625.       }
  626.       sprintf(sz, " (%.2f%c %.2f%c)",
  627.         yV, yA < 0.0 ? 's' : 'n', xV, xA > 0.0 ? 'e' : 'w'); PrintSz(sz);
  628.     } else
  629.       PrintAltitude(azialt[i]);
  630.     PrintL();
  631.   }
  632.   if (occurcount == 0)
  633.     PrintSz("No horizon events found.\n");
  634.  
  635.   /* Recompute original chart placements as we've overwritten them. */
  636.  
  637.   ciCore = ciMain; ciTwin = ciT;
  638.   us.fSidereal = fT;
  639.   CastChart(fTrue);
  640. }
  641.  
  642.  
  643. /* Print out an ephemeris - the positions of the planets (at the time in the */
  644. /* current chart) each day during a specified month, as done with the -E     */
  645. /* switch. Display the ephemeris for the whole year if -Ey is in effect.     */
  646.  
  647. void ChartEphemeris()
  648. {
  649.   char sz[cchSzDef];
  650.   int yea, yea1, yea2, mon, mon1, mon2, daysiz, i, j, s, d, m;
  651.  
  652.   /* If -Ey is in effect, then loop through all months in the whole year. */
  653.  
  654.   if (us.nEphemYears) {
  655.     yea1 = Yea; yea2 = yea1 + us.nEphemYears - 1; mon1 = 1; mon2 = 12;
  656.   } else {
  657.     yea1 = yea2 = Yea; mon1 = mon2 = Mon;
  658.   }
  659.  
  660.   /* Loop through the year or years in question. */
  661.  
  662.   for (yea = yea1; yea <= yea2; yea++)
  663.  
  664.   /* Loop through the month or months in question, printing each ephemeris. */
  665.  
  666.   for (mon = mon1; mon <= mon2; mon++) {
  667.     daysiz = DayInMonth(mon, yea);
  668.     PrintSz(us.fEuroDate ? "Dy/Mo/Yr" : "Mo/Dy/Yr");
  669.     for (j = 1; j <= cObj; j++) {
  670.       if (!ignore[j] && FThing(j)) {
  671.         sprintf(sz, "  %s%c%c%c%c", is.fSeconds ? "  " : "", chObj3(j),
  672.           szObjName[j][3] != 0 ? szObjName[j][3] : ' '); PrintSz(sz);
  673.         PrintTab(' ', us.fParallel ? 2 + is.fSeconds : 1 + 3*is.fSeconds);
  674.       }
  675.     }
  676.     PrintL();
  677.     for (i = 1; i <= daysiz; i = AddDay(mon, i, yea, 1)) {
  678.  
  679.       /* Loop through each day in the month, casting a chart for that day. */
  680.  
  681.       SetCI(ciCore, mon, i, yea, Tim, Dst, Zon, Lon, Lat);
  682.       CastChart(fTrue);
  683.       PrintSz(SzDate(mon, i, yea, -1));
  684.       PrintCh(' ');
  685.       for (j = 0; j <= cObj; j++)
  686.         if (!FIgnore(j) && FThing(j)) {
  687.           if (!us.fParallel) {
  688.             if (is.fSeconds)
  689.               PrintZodiac(planet[j]);
  690.             else {
  691.               AnsiColor(kObjA[j]);
  692.               s = SFromZ(planet[j]);
  693.               d = (int)planet[j] - (s-1)*30;
  694.               m = (int)(RFract(planet[j])*60.0);
  695.               sprintf(sz, "%2d%s%02d", d, szSignAbbrev[s], m); PrintSz(sz);
  696.             }
  697.           } else {
  698.             AnsiColor(kObjA[j]);
  699.             PrintAltitude(planetalt[j]);
  700.           }
  701.           PrintCh((char)(ret[j] >= 0.0 ? ' ' : '.'));
  702.         }
  703.       PrintL();
  704.       AnsiColor(kDefault);
  705.     }
  706.     if (mon < mon2 || yea < yea2)
  707.       PrintL();
  708.   }
  709.  
  710.   ciCore = ciMain;    /* Recast original chart. */
  711.   CastChart(fTrue);
  712. }
  713.  
  714. /* charts3.c */
  715.